home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / CRowVector.cc < prev    next >
C/C++ Source or Header  |  1997-03-07  |  16KB  |  803 lines
  1. // RowVector manipulations.
  2. /*
  3.  
  4. Copyright (C) 1996 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  21.  
  22. */
  23.  
  24. #if defined (__GNUG__)
  25. #pragma implementation
  26. #endif
  27.  
  28. #ifdef HAVE_CONFIG_H
  29. #include <config.h>
  30. #endif
  31.  
  32. #include <iostream.h>
  33.  
  34. #include "f77-fcn.h"
  35. #include "lo-error.h"
  36. #include "mx-base.h"
  37. #include "mx-inlines.cc"
  38. #include "oct-cmplx.h"
  39.  
  40. // Fortran functions we call.
  41.  
  42. extern "C"
  43. {
  44.   int F77_FCN (zgemv, ZGEMV) (const char*, const int&, const int&,
  45.                   const Complex&, const Complex*,
  46.                   const int&, const Complex*, const int&,
  47.                   const Complex&, Complex*, const int&,
  48.                   long);
  49. }
  50.  
  51. // Complex Row Vector class
  52.  
  53. ComplexRowVector::ComplexRowVector (const RowVector& a)
  54.   : MArray<Complex> (a.length ())
  55. {
  56.   for (int i = 0; i < length (); i++)
  57.     elem (i) = a.elem (i);
  58. }
  59.  
  60. bool
  61. ComplexRowVector::operator == (const ComplexRowVector& a) const
  62. {
  63.   int len = length ();
  64.   if (len != a.length ())
  65.     return 0;
  66.   return equal (data (), a.data (), len);
  67. }
  68.  
  69. bool
  70. ComplexRowVector::operator != (const ComplexRowVector& a) const
  71. {
  72.   return !(*this == a);
  73. }
  74.  
  75. // destructive insert/delete/reorder operations
  76.  
  77. ComplexRowVector&
  78. ComplexRowVector::insert (const RowVector& a, int c)
  79. {
  80.   int a_len = a.length ();
  81.   if (c < 0 || c + a_len > length ())
  82.     {
  83.       (*current_liboctave_error_handler) ("range error for insert");
  84.       return *this;
  85.     }
  86.  
  87.   for (int i = 0; i < a_len; i++)
  88.     elem (c+i) = a.elem (i);
  89.  
  90.   return *this;
  91. }
  92.  
  93. ComplexRowVector&
  94. ComplexRowVector::insert (const ComplexRowVector& a, int c)
  95. {
  96.   int a_len = a.length ();
  97.   if (c < 0 || c + a_len > length ())
  98.     {
  99.       (*current_liboctave_error_handler) ("range error for insert");
  100.       return *this;
  101.     }
  102.  
  103.   for (int i = 0; i < a_len; i++)
  104.     elem (c+i) = a.elem (i);
  105.  
  106.   return *this;
  107. }
  108.  
  109. ComplexRowVector&
  110. ComplexRowVector::fill (double val)
  111. {
  112.   int len = length ();
  113.   if (len > 0)
  114.     for (int i = 0; i < len; i++)
  115.       elem (i) = val;
  116.   return *this;
  117. }
  118.  
  119. ComplexRowVector&
  120. ComplexRowVector::fill (const Complex& val)
  121. {
  122.   int len = length ();
  123.   if (len > 0)
  124.     for (int i = 0; i < len; i++)
  125.       elem (i) = val;
  126.   return *this;
  127. }
  128.  
  129. ComplexRowVector&
  130. ComplexRowVector::fill (double val, int c1, int c2)
  131. {
  132.   int len = length ();
  133.   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
  134.     {
  135.       (*current_liboctave_error_handler) ("range error for fill");
  136.       return *this;
  137.     }
  138.  
  139.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  140.  
  141.   for (int i = c1; i <= c2; i++)
  142.     elem (i) = val;
  143.  
  144.   return *this;
  145. }
  146.  
  147. ComplexRowVector&
  148. ComplexRowVector::fill (const Complex& val, int c1, int c2)
  149. {
  150.   int len = length ();
  151.   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
  152.     {
  153.       (*current_liboctave_error_handler) ("range error for fill");
  154.       return *this;
  155.     }
  156.  
  157.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  158.  
  159.   for (int i = c1; i <= c2; i++)
  160.     elem (i) = val;
  161.  
  162.   return *this;
  163. }
  164.  
  165. ComplexRowVector
  166. ComplexRowVector::append (const RowVector& a) const
  167. {
  168.   int len = length ();
  169.   int nc_insert = len;
  170.   ComplexRowVector retval (len + a.length ());
  171.   retval.insert (*this, 0);
  172.   retval.insert (a, nc_insert);
  173.   return retval;
  174. }
  175.  
  176. ComplexRowVector
  177. ComplexRowVector::append (const ComplexRowVector& a) const
  178. {
  179.   int len = length ();
  180.   int nc_insert = len;
  181.   ComplexRowVector retval (len + a.length ());
  182.   retval.insert (*this, 0);
  183.   retval.insert (a, nc_insert);
  184.   return retval;
  185. }
  186.  
  187. ComplexColumnVector
  188. ComplexRowVector::hermitian (void) const
  189. {
  190.   int len = length ();
  191.   return ComplexColumnVector (conj_dup (data (), len), len);
  192. }
  193.  
  194. ComplexColumnVector
  195. ComplexRowVector::transpose (void) const
  196. {
  197.   return ComplexColumnVector (*this);
  198. }
  199.  
  200. ComplexRowVector
  201. conj (const ComplexRowVector& a)
  202. {
  203.   int a_len = a.length ();
  204.   ComplexRowVector retval;
  205.   if (a_len > 0)
  206.     retval = ComplexRowVector (conj_dup (a.data (), a_len), a_len);
  207.   return retval;
  208. }
  209.  
  210. // resize is the destructive equivalent for this one
  211.  
  212. ComplexRowVector
  213. ComplexRowVector::extract (int c1, int c2) const
  214. {
  215.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  216.  
  217.   int new_c = c2 - c1 + 1;
  218.  
  219.   ComplexRowVector result (new_c);
  220.  
  221.   for (int i = 0; i < new_c; i++)
  222.     result.elem (i) = elem (c1+i);
  223.  
  224.   return result;
  225. }
  226.  
  227. // row vector by row vector -> row vector operations
  228.  
  229. ComplexRowVector&
  230. ComplexRowVector::operator += (const RowVector& a)
  231. {
  232.   int len = length ();
  233.  
  234.   int a_len = a.length ();
  235.  
  236.   if (len != a_len)
  237.     {
  238.       gripe_nonconformant ("operator +=", len, a_len);
  239.       return *this;
  240.     }
  241.  
  242.   if (len == 0)
  243.     return *this;
  244.  
  245.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  246.  
  247.   add2 (d, a.data (), len);
  248.   return *this;
  249. }
  250.  
  251. ComplexRowVector&
  252. ComplexRowVector::operator -= (const RowVector& a)
  253. {
  254.   int len = length ();
  255.  
  256.   int a_len = a.length ();
  257.  
  258.   if (len != a_len)
  259.     {
  260.       gripe_nonconformant ("operator -=", len, a_len);
  261.       return *this;
  262.     }
  263.  
  264.   if (len == 0)
  265.     return *this;
  266.  
  267.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  268.  
  269.   subtract2 (d, a.data (), len);
  270.   return *this;
  271. }
  272.  
  273. ComplexRowVector&
  274. ComplexRowVector::operator += (const ComplexRowVector& a)
  275. {
  276.   int len = length ();
  277.  
  278.   int a_len = a.length ();
  279.  
  280.   if (len != a_len)
  281.     {
  282.       gripe_nonconformant ("operator +=", len, a_len);
  283.       return *this;
  284.     }
  285.  
  286.   if (len == 0)
  287.     return *this;
  288.  
  289.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  290.  
  291.   add2 (d, a.data (), len);
  292.   return *this;
  293. }
  294.  
  295. ComplexRowVector&
  296. ComplexRowVector::operator -= (const ComplexRowVector& a)
  297. {
  298.   int len = length ();
  299.  
  300.   int a_len = a.length ();
  301.  
  302.   if (len != a_len)
  303.     {
  304.       gripe_nonconformant ("operator -=", len, a_len);
  305.       return *this;
  306.     }
  307.  
  308.   if (len == 0)
  309.     return *this;
  310.  
  311.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  312.  
  313.   subtract2 (d, a.data (), len);
  314.   return *this;
  315. }
  316.  
  317. // row vector by scalar -> row vector operations
  318.  
  319. ComplexRowVector
  320. operator + (const ComplexRowVector& v, double s)
  321. {
  322.   int len = v.length ();
  323.   return ComplexRowVector (add (v.data (), len, s), len);
  324. }
  325.  
  326. ComplexRowVector
  327. operator - (const ComplexRowVector& v, double s)
  328. {
  329.   int len = v.length ();
  330.   return ComplexRowVector (subtract (v.data (), len, s), len);
  331. }
  332.  
  333. ComplexRowVector
  334. operator * (const ComplexRowVector& v, double s)
  335. {
  336.   int len = v.length ();
  337.   return ComplexRowVector (multiply (v.data (), len, s), len);
  338. }
  339.  
  340. ComplexRowVector
  341. operator / (const ComplexRowVector& v, double s)
  342. {
  343.   int len = v.length ();
  344.   return ComplexRowVector (divide (v.data (), len, s), len);
  345. }
  346.  
  347. ComplexRowVector
  348. operator + (const RowVector& v, const Complex& s)
  349. {
  350.   int len = v.length ();
  351.   return ComplexRowVector (add (v.data (), len, s), len);
  352. }
  353.  
  354. ComplexRowVector
  355. operator - (const RowVector& v, const Complex& s)
  356. {
  357.   int len = v.length ();
  358.   return ComplexRowVector (subtract (v.data (), len, s), len);
  359. }
  360.  
  361. ComplexRowVector
  362. operator * (const RowVector& v, const Complex& s)
  363. {
  364.   int len = v.length ();
  365.   return ComplexRowVector (multiply (v.data (), len, s), len);
  366. }
  367.  
  368. ComplexRowVector
  369. operator / (const RowVector& v, const Complex& s)
  370. {
  371.   int len = v.length ();
  372.   return ComplexRowVector (divide (v.data (), len, s), len);
  373. }
  374.  
  375. // scalar by row vector -> row vector operations
  376.  
  377. ComplexRowVector
  378. operator + (double s, const ComplexRowVector& a)
  379. {
  380.   int a_len = a.length ();
  381.   return ComplexRowVector (add (a.data (), a_len, s), a_len);
  382. }
  383.  
  384. ComplexRowVector
  385. operator - (double s, const ComplexRowVector& a)
  386. {
  387.   int a_len = a.length ();
  388.   return ComplexRowVector (subtract (s, a.data (), a_len), a_len);
  389. }
  390.  
  391. ComplexRowVector
  392. operator * (double s, const ComplexRowVector& a)
  393. {
  394.   int a_len = a.length ();
  395.   return ComplexRowVector (multiply (a.data (), a_len, s), a_len);
  396. }
  397.  
  398. ComplexRowVector
  399. operator / (double s, const ComplexRowVector& a)
  400. {
  401.   int a_len = a.length ();
  402.   return ComplexRowVector (divide (s, a.data (), a_len), a_len);
  403. }
  404.  
  405. ComplexRowVector
  406. operator + (const Complex& s, const RowVector& a)
  407. {
  408.   int a_len = a.length ();
  409.   return ComplexRowVector (add (a.data (), a_len, s), a_len);
  410. }
  411.  
  412. ComplexRowVector
  413. operator - (const Complex& s, const RowVector& a)
  414. {
  415.   int a_len = a.length ();
  416.   return ComplexRowVector (subtract (s, a.data (), a_len), a_len);
  417. }
  418.  
  419. ComplexRowVector
  420. operator * (const Complex& s, const RowVector& a)
  421. {
  422.   int a_len = a.length ();
  423.   return ComplexRowVector (multiply (a.data (), a_len, s), a_len);
  424. }
  425.  
  426. ComplexRowVector
  427. operator / (const Complex& s, const RowVector& a)
  428. {
  429.   int a_len = a.length ();
  430.   return ComplexRowVector (divide (s, a.data (), a_len), a_len);
  431. }
  432.  
  433. // row vector by matrix -> row vector
  434.  
  435. ComplexRowVector
  436. operator * (const ComplexRowVector& v, const ComplexMatrix& a)
  437. {
  438.   ComplexRowVector retval;
  439.  
  440.   int len = v.length ();
  441.  
  442.   int a_nr = a.rows ();
  443.   int a_nc = a.cols ();
  444.  
  445.   if (a_nr != len)
  446.     gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
  447.   else
  448.     {
  449.       int a_nr = a.rows ();
  450.       int a_nc = a.cols ();
  451.  
  452.       if (len == 0)
  453.     retval.resize (a_nc, 0.0);
  454.       else
  455.     {
  456.       // Transpose A to form A'*x == (x'*A)'
  457.  
  458.       int ld = a_nr;
  459.  
  460.       retval.resize (a_nc);
  461.       Complex *y = retval.fortran_vec ();
  462.  
  463.       F77_XFCN (zgemv, ZGEMV, ("T", a_nr, a_nc, 1.0, a.data (),
  464.                    ld, v.data (), 1, 0.0, y, 1, 1L));
  465.  
  466.       if (f77_exception_encountered)
  467.         (*current_liboctave_error_handler)
  468.           ("unrecoverable error in zgemv");
  469.     }
  470.     }
  471.  
  472.   return retval;
  473. }
  474.  
  475. ComplexRowVector
  476. operator * (const RowVector& v, const ComplexMatrix& a)
  477. {
  478.   ComplexRowVector tmp (v);
  479.   return tmp * a;
  480. }
  481.  
  482. // row vector by row vector -> row vector operations
  483.  
  484. ComplexRowVector
  485. operator + (const ComplexRowVector& v, const RowVector& a)
  486. {
  487.   int len = v.length ();
  488.  
  489.   int a_len = a.length ();
  490.  
  491.   if (len != a_len)
  492.     {
  493.       gripe_nonconformant ("operator +", len, a_len);
  494.       return ComplexRowVector ();
  495.     }
  496.  
  497.   if (len == 0)
  498.     return ComplexRowVector (0);
  499.  
  500.   return ComplexRowVector (add (v.data (), a.data (), len), len);
  501. }
  502.  
  503. ComplexRowVector
  504. operator - (const ComplexRowVector& v, const RowVector& a)
  505. {
  506.   int len = v.length ();
  507.  
  508.   int a_len = a.length ();
  509.  
  510.   if (len != a_len)
  511.     {
  512.       gripe_nonconformant ("operator -", len, a_len);
  513.       return ComplexRowVector ();
  514.     }
  515.  
  516.   if (len == 0)
  517.     return ComplexRowVector (0);
  518.  
  519.   return ComplexRowVector (subtract (v.data (), a.data (), len), len);
  520. }
  521.  
  522. ComplexRowVector
  523. operator + (const RowVector& v, const ComplexRowVector& a)
  524. {
  525.   int len = v.length ();
  526.  
  527.   int a_len = a.length ();
  528.  
  529.   if (len != a_len)
  530.     {
  531.       gripe_nonconformant ("operator +", len, a_len);
  532.       return ComplexRowVector ();
  533.     }
  534.  
  535.   if (len == 0)
  536.     return ComplexRowVector (0);
  537.  
  538.   return ComplexRowVector (add (v.data (), a.data (), len), len);
  539. }
  540.  
  541. ComplexRowVector
  542. operator - (const RowVector& v, const ComplexRowVector& a)
  543. {
  544.   int len = v.length ();
  545.  
  546.   int a_len = a.length ();
  547.  
  548.   if (len != a_len)
  549.     {
  550.       gripe_nonconformant ("operator -", len, a_len);
  551.       return ComplexRowVector ();
  552.     }
  553.  
  554.   if (len == 0)
  555.     return ComplexRowVector (0);
  556.  
  557.   return ComplexRowVector (subtract (v.data (), a.data (), len), len);
  558. }
  559.  
  560. ComplexRowVector
  561. product (const ComplexRowVector& v, const RowVector& a)
  562. {
  563.   int len = v.length ();
  564.  
  565.   int a_len = a.length ();
  566.  
  567.   if (len != a_len)
  568.     {
  569.       gripe_nonconformant ("product", len, a_len);
  570.       return ComplexRowVector ();
  571.     }
  572.  
  573.   if (len == 0)
  574.     return ComplexRowVector (0);
  575.  
  576.   return ComplexRowVector (multiply (v.data (), a.data (), len), len);
  577. }
  578.  
  579. ComplexRowVector
  580. quotient (const ComplexRowVector& v, const RowVector& a)
  581. {
  582.   int len = v.length ();
  583.  
  584.   int a_len = a.length ();
  585.  
  586.   if (len != a_len)
  587.     {
  588.       gripe_nonconformant ("quotient", len, a_len);
  589.       return ComplexRowVector ();
  590.     }
  591.  
  592.   if (len == 0)
  593.     return ComplexRowVector (0);
  594.  
  595.   return ComplexRowVector (divide (v.data (), a.data (), len), len);
  596. }
  597.  
  598. ComplexRowVector
  599. product (const RowVector& v, const ComplexRowVector& a)
  600. {
  601.   int len = v.length ();
  602.  
  603.   int a_len = a.length ();
  604.  
  605.   if (len != a_len)
  606.     {
  607.       gripe_nonconformant ("product", len, a_len);
  608.       return ComplexRowVector ();
  609.     }
  610.  
  611.   if (len == 0)
  612.     return ComplexRowVector (0);
  613.  
  614.   return ComplexRowVector (multiply (v.data (), a.data (), len), len);
  615. }
  616.  
  617. ComplexRowVector
  618. quotient (const RowVector& v, const ComplexRowVector& a)
  619. {
  620.   int len = v.length ();
  621.  
  622.   int a_len = a.length ();
  623.  
  624.   if (len != a_len)
  625.     {
  626.       gripe_nonconformant ("quotient", len, a_len);
  627.       return ComplexRowVector ();
  628.     }
  629.  
  630.   if (len == 0)
  631.     return ComplexRowVector (0);
  632.  
  633.   return ComplexRowVector (divide (v.data (), a.data (), len), len);
  634. }
  635.  
  636. // other operations
  637.  
  638. ComplexRowVector
  639. ComplexRowVector::map (c_c_Mapper f) const
  640. {
  641.   ComplexRowVector b (*this);
  642.   return b.apply (f);
  643. }
  644.  
  645. RowVector
  646. ComplexRowVector::map (d_c_Mapper f) const
  647. {
  648.   const Complex *d = data ();
  649.  
  650.   int len = length ();
  651.  
  652.   RowVector retval (len);
  653.  
  654.   double *r = retval.fortran_vec ();
  655.  
  656.   for (int i = 0; i < len; i++)
  657.     r[i] = f (d[i]);
  658.  
  659.   return retval;
  660. }
  661.  
  662. ComplexRowVector&
  663. ComplexRowVector::apply (c_c_Mapper f)
  664. {
  665.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  666.  
  667.   for (int i = 0; i < length (); i++)
  668.     d[i] = f (d[i]);
  669.  
  670.   return *this;
  671. }
  672.  
  673. Complex
  674. ComplexRowVector::min (void) const
  675. {
  676.   int len = length ();
  677.   if (len == 0)
  678.     return Complex (0.0);
  679.  
  680.   Complex res = elem (0);
  681.   double absres = abs (res);
  682.  
  683.   for (int i = 1; i < len; i++)
  684.     if (abs (elem (i)) < absres)
  685.       {
  686.     res = elem (i);
  687.     absres = abs (res);
  688.       }
  689.  
  690.   return res;
  691. }
  692.  
  693. Complex
  694. ComplexRowVector::max (void) const
  695. {
  696.   int len = length ();
  697.   if (len == 0)
  698.     return Complex (0.0);
  699.  
  700.   Complex res = elem (0);
  701.   double absres = abs (res);
  702.  
  703.   for (int i = 1; i < len; i++)
  704.     if (abs (elem (i)) > absres)
  705.       {
  706.     res = elem (i);
  707.     absres = abs (res);
  708.       }
  709.  
  710.   return res;
  711. }
  712.  
  713. // i/o
  714.  
  715. ostream&
  716. operator << (ostream& os, const ComplexRowVector& a)
  717. {
  718. //  int field_width = os.precision () + 7;
  719.   for (int i = 0; i < a.length (); i++)
  720.     os << " " /* setw (field_width) */ << a.elem (i);
  721.   return os;
  722. }
  723.  
  724. istream&
  725. operator >> (istream& is, ComplexRowVector& a)
  726. {
  727.   int len = a.length();
  728.  
  729.   if (len < 1)
  730.     is.clear (ios::badbit);
  731.   else
  732.     {
  733.       Complex tmp;
  734.       for (int i = 0; i < len; i++)
  735.         {
  736.           is >> tmp;
  737.           if (is)
  738.             a.elem (i) = tmp;
  739.           else
  740.             break;
  741.         }
  742.     }
  743.   return is;
  744. }
  745.  
  746. // row vector by column vector -> scalar
  747.  
  748. // row vector by column vector -> scalar
  749.  
  750. Complex
  751. operator * (const ComplexRowVector& v, const ColumnVector& a)
  752. {
  753.   ComplexColumnVector tmp (a);
  754.   return v * tmp;
  755. }
  756.  
  757. Complex
  758. operator * (const ComplexRowVector& v, const ComplexColumnVector& a)
  759. {
  760.   int len = v.length ();
  761.  
  762.   int a_len = a.length ();
  763.  
  764.   if (len != a_len)
  765.     {
  766.       gripe_nonconformant ("operator *", len, a_len);
  767.       return 0.0;
  768.     }
  769.  
  770.   Complex retval (0.0, 0.0);
  771.  
  772.   for (int i = 0; i < len; i++)
  773.     retval += v.elem (i) * a.elem (i);
  774.  
  775.   return retval;
  776. }
  777.  
  778. // other operations
  779.  
  780. ComplexRowVector
  781. linspace (const Complex& x1, const Complex& x2, int n)
  782. {
  783.   ComplexRowVector retval;
  784.  
  785.   if (n > 0)
  786.     {
  787.       retval.resize (n);
  788.       Complex delta = (x2 - x1) / (n - 1);
  789.       retval.elem (0) = x1;
  790.       for (int i = 1; i < n-1; i++)
  791.     retval.elem (i) = x1 + i * delta;
  792.       retval.elem (n-1) = x2;
  793.     }
  794.  
  795.   return retval;
  796. }
  797.  
  798. /*
  799. ;;; Local Variables: ***
  800. ;;; mode: C++ ***
  801. ;;; End: ***
  802. */
  803.